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We have compared the results of our approximation scheme, the composite operator 
method, for the double occupancy and the internal energy of the two-dimensional 
Hubbard model with numerical data obtained by means of the Lanczos^ and quan- 
tum Monte Carlo schemes^. The agreement is very good at both half- filling and 
away from it showing how reliable is the approximation scheme. 

1 Introduction 

The relevance of the Hubbard model^ as minimal model for describing the high- 
Tc cuprate superconductors'* and, generally, many strongly correlated materials^ is 
very well established. The model is thought to have a very rich phenomenology: 
a variety of spin and charge orderings^ (ferro- and antiferro- magnetic included), 
Mott transitions (of the Slater, Heisenberg and Hubbard types) driven by both 
filling and coupling^' ^, non- Fermi liquid dynamics, superconductivity caused by: 
magnons, static and dynamic charge ordering, proximity to quantum critical points 
of different origin. In particular, strong antiferromagnetic correlations are present 
at half-filling (and for low doping) and low temperatures^. 

According to this, it is very relevant to have reliable solutions of this model for 
the variety of boundary conditions under which it can be studied. Actually, we 
know the exact solution only in one dimension thanks to the Bethe Ansatz^. Some 
other few exact results are known in quite special limits, but the more relevant 
questions are still unsolved as the exact solution in two dimensions (which seems 
to be the case for the majority of the emergent highly interacting materials) is 
still missing. Owing to this, a huge number of approximation schemes have been 
proposed since the very beginning (Hubbard presented the model together with an 
approximate solution known today as Hubbard I) and their reliability is still under 
test as many physical interpretations are based on the results obtained by them. 

In the last decades we have seen the birth and development of many numerical 
methods (exact diagonalization, Lanczos, quantum Monte Carlo, ...) for studying 
finite clusters of bigger and bigger size. The results of these numerical methods 
are extremely important for the development of analytical schemes as they provide 
the possibility to execute unbiased tests. The numerical data can be interpreted as 
the experimental results relative to a specific model and any comparison with an 
analytical approximation scheme can be directly contrasted with its reliability. On 
the contrary, comparisons with real experimental data suffer of an high degree of 
uncertainty regarding the capability of the chosen model to capture the essential 
physical of the material under analysis. 

In this last decade, we have been developing an approximation scheme, namely 
the composite operator method^' which proved to be quite powerful to de- 
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scribe the physics of many correlated models and the properties of some materi- 
als. The main peculiarities of this method are: the use of symmetry constraints 
(Pauli principle and Ward-Takahashi identities) to fix the representation where the 
Green's functions are realized, the freedom in choosing the composite fields for the 
operatorial basis around which we wish to construct perturbative solutions. The 
algebra imposes many relations between operators which have to be fulfilled also 
at the level of averages. According to this, it is possible to individuate a set of 
self-consistent equations which allows to fix the unknown parameters coming in the 
calculations and, consequently, to instruct the solution regarding the Hilbert space 
in which the composite operators live (i.e., the representation). In particular, we 
see the unknown parameters as a necessity and not as an accident as many other 
approximation schemes do. In the last years, we have acquired much experience in 
choosing the composite operators more effective in describing the physical property 
of a system and we have collected a certain number of recipes which can be used 
for specific class of models. 

In this manuscript, we wish to present some new comparisons between numerical 
results and composite operator method ones for two local quantities: the double 
occupancy and the internal energy. 



2 Results 

The Hubbard model is described by the following Hamiltonian 

H = J^{Ui-f, <5ij) ct {i) c (i) + C/ ^ nt (i) (i) (1) 

ij i 

where (i) = (i) , c| {ifj is the electronic creation operator in spinorial notation 

at the site i [i = (i, t)] and (i) = (i) c„ (i) is the number operator for spin a 
at the site i; /x is the chemical potential and U is the on-site Coulomb repulsion. 

The matrix describes the nearest-neighbor hopping; in the 2D case we have 
tij = — 4taij, where 



l^e"'('-J)a(k) (2) 



N 



k 



is the projector on the nearest-neighbor sites and a (k) = ^ [cos {k^ a) + cos {ky a)] 

and a is the lattice parameter. 

We choose the following fermionic basis 

(3) 



T]{i) 

where ^{i) = [1 — n {i)] c {i) and r]{i) = n {i) c (i) are the Hubbard operators. ^ (i) 
satisfies the following equation of motion 
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where c" (i, t) = X]j '^ij 0' ^) ^^"^ ''^ (*) = ^cr'^ (i) c" (?) + c (i) [c^'" (i) c (i)] . 
nju(i) = c^{i) (J^ c{i) are the charge (/i = 0) and spin (y^ = 1, 2, 3 ) density operators, 
with (jfj, = (1,1?), = {—1,5) and a are the Pauh matrices. 
Let us project the source J (i) on the chosen basis 

J(i,t)-5^£(i,j)*(j,i) (5) 
j 

Accordingly, the energy matrix e (i, j) is defined through the equation 

m(i,j) = ^£(i,l)/(l,j) (6) 
1 

where the m-matrix and the normalization matrix I have the following definitions 

m(i,j) = ({J(M),vl/t(j,t)}) (7) 
7(i,j) = <{*(i,i),*t(j,t)}) (8) 

After Eq. (5), the Fourier transform of the thermal single-particle retarded Green's 
function G = {R {i) (j)] ) satisfies the following equation 

[a;-e(k)]G(k,a;)=7(k) (9) 

The straightforward application of this scheme*' ^""^'^ gives that, in the para- 
magnetic phase, /(k) has diagonal form with In = 1— n/2 and I22 = 'n/2 
{{ric, (i)) = f ) and that the m-matrix depends on three parameters: the chemi- 
cal potential /U and two correlators 

A={^"ii)^Hi))-{'n'^ii)vHi)) (10) 

P=lK (^) (0) - ([CT (i) Ci (iTcl (i) c\ (i)) (11) 

The three parameters /x, A and p are functions of the external parameters n, T and 
U and can be determined self-consistently through a set of three coupled equations 

n = 2[l-(c(i)ct(i))] 

A = (e" (0 (i)) - (i) (i)) (12) 

iCiz) ,?t (,))=o 

The first equation fixes the particle number, the second one comes from the defi- 
nition of A and the third one assures that the solution respects the Pauli principle 
(i.e., the local algebra)'^''. 

In Fig. 1 we report the behavior of the internal energy per site Eg = {H) /N as 
a function of the Coulomb interaction U at half- filling and n = 8/9. The Lanczos 
numerical data have been taken from Rcf. ^. The agreement at half- filling is very 
good. We have got an exceptionally good agreement at half-filling for the one- 
dimensional Hubbard model too^^. In that case, we have compared our results 
with the exact Bethe ansatz solution. Away from half-filling, the agreement is 
less satisfactory, but, as we will see from the double occupancy comparison, the 
difference is quite constant as function of U and can be due to finite size effects, as 
surely is the difference at U = 0. 
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Figure 1 . Internal energy per site Eg as a function of the Coulomb interaction U for T = and 
n = 1 and 8/9. The Lanczos data are taken from Ref. ^. 
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Figure 2. Double occupancy D as a function of the Coulomb interaction U for T = 1/16 and 
n = 1. The quantum Monte Carlo data are taJsen from Ref. ^. 

The double oecupaney D = {npi i), whieh ean be computed by means of ther- 
modynamics through the following formula D = ^ {F is the free energy; for T = 
we have D = ) , is reported as a function of the Coulomb interaction U in Fig. 2 
for T = 1/16 and n = 1 and in Fig. 3 for T = and n = 8/9. The Lanczos numer- 
ical data have again been taken from Ref. ^; the quantum Monte Carlo ones from 
Ref. ^. The agreement is very good for both values of filling. For half- filling, we 

have also reported the limiting behavior of our solution: D where d is the 

spatial dimension of the system^^. The numerical data seem to follow this behavior 
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Figure 3. Double occupancy D as a function of the Coulomb interaction U for T = and n = 8/9. 
The Lanczos data are taJsen from Ref. ^. 

also for not so large value of the Coulomb repulsfon. This is probably due to finite 

size effects. For n = 8/9, the very good agreement supports our interpretation of 
the comparison for the corresponding energy. 

3 Conclusions 

The positive comparisons with the Bethe ansatz exact solution for the one- 
dimensional case and with the numerical results for the two-dimensional one show 
how reliable is our approximation scheme to study the Hubbard model. In particu- 
lar, the agreement reported in this manuscript with the Lanczos and the quantum 
Monte Carlo data at both half-filling and away from it is very good and shows our 
capability to properly describe both the charge and the spin fluctuation scales of 
energy. 
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